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Abstract. The mesochronic velocity is the average of the velocity field along 
trajectories generated by the same velocity field over a time interval of finite 
duration. In this paper we classify initial conditions of trajectories evolving 
in incompressible vector fields according to the character of motion of mate¬ 
rial around the trajectory. In particular, we provide calculations that can be 
used to determine the number of expanding directions and the presence of 
rotation from the characteristic polynomial of the Jacobian matrix of meso¬ 
chronic velocity. In doing so, we show that (a) the mesochronic velocity can 
be used to characterize dynamical deformation of three-dimensional volumes, 
(b) the resulting mesochronic analysis is a finite-time extension of the Okubo- 
Weiss-Chong analysis of incompressible velocity fields, (c) the two-dimensional 
mesochronic analysis from Mezic et al. “A New Mixing Diagnostic and Gulf Oil 
Spill Movement”, Science 330, (2010), 486—489, extends to three-dimensional 
state spaces. Theoretical considerations are further supported by numerical 
computations performed for a dynamical system arising in fluid mechanics, the 
unsteady Arnold-Beltrami-Childress (ABC) flow. 
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1. Introduction 

Chaotic advection is the theory of material transport in fluids based in dynamical 
systems theory. [3, 52,65] It is largely rooted in analysis of geometric structures in 
flows with a simple time-dependence, time-autonomous or periodic. Since the real¬ 
ization that even flows with a complicated time dependence, e.g., turbulent flows, 
possess organized Lagrangian structures, it has become increasingly important to 
detect geometric structures analogous to invariant manifolds in steady flows. In 
particular, detection of structures using only finite-time information about the flow 
has been seen as the most practically-useful direction. [2,21,28,32,56,59,60] 

The need for analysis of geometric structures that organize advection is not 
purely academic. Transport of material by fluid flows played a crucial role in the 
fallout from several recent catastrophic events, namely, the volcanic eruption of 
Eyjafjallajokull (2010), the Deepwater Horizon oil spill (2010), and the nuclear 
disaster in Fukushima (2011). These events highlighted how important it is to detect 
organizing geometric structures in (near) real time from data, either measured or 
generated by detailed simulation models; consequently, such problems have become 
a very active intersection of dynamical systems and fluid dynamics. 

The problem of identifying geometric structures in flows has resulted in several 
approaches, each focusing on a somewhat-different structure as the objective of its 
analysis. 

The theory of Lagrangian Coherent Struetures ]33] (LCS) identifies barriers that 
organize the transport in flows with complex time-dependence. Initially, LCS were 
closely associated with computation of Finite-Time Lyapunov Exponent fields ]60]; 
more recently, they have been re-formulated using a variational principle [5,29,31], 
which defines them as certain geodesic lines of the local deformation held induced 
by the fluid flow. This new definition allows a finer classification of LCS, both in two 
and three dimensions, based on the type of deformation, e.g., hyperbolic, elliptic, 
corresponding to different behaviors of fluid parcels in the flow. The recent review 
by Haller [30] gives a detailed coverage of the current state-of-the-art of techniques 
centered around LCS. 

LCS and associated theories mostly focused on magnitude of non-rotational de¬ 
formation in the flow. The rotation deformation has been classically studied by 
Poincare in topological dynamics and Ruelle [58] in ergodic theory; recently, a 
Finite-Time Rotation Number [61] has been proposed as a useful quantity in anal¬ 
ysis of flows. At the closing of this manuscript we were also notified of the recent 
work by Farazmand and Haller [16], working along the same lines. 

In the effort to study the “stretch-and-fold” mechanism for chaos in finite time, 
most studies focus on the first-order “stretch”. Folding in finite time has not re¬ 
ceived direct attention; an exception is the study of the Finite-Time Curvature 
Fields [37-39]. At this point, structures observed through all these methods have 
been connected mostly on phenomenological basis, through comparison of visual¬ 
izations, showing considerable overlap between observed structures but, also, non- 
negligible differences. 

Magnitudes of the local material deformation are typically estimated by pro¬ 
cessing velocity gradients; they cannot be precisely computed in the absence of the 
detailed data about the velocity held, e.g., when the system is sampled by sparse tra¬ 
jectories only. In sparsely-sampled planar systems, trajectories can be represented 
by space-time braids — extremely-reduced, symbolic representations of trajectories. 
The resulting approach, known as braid dynamics [2, 6,62] has been successful in 
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providing lower bounds on the amount of topological deformation present in the 
flow, in limited-data settings. The obtained bounds have been used both in de¬ 
sign and analysis of the material advection; unfortunately, there are currently no 
extensions of braid dynamics to three-dimensional flows. 

Instead of looking for barriers to transport, as is the case with the LCS theory, 
the theory of almost-invariant sets identifies a collection of sets, fixed in space, such 
that the material placed in them does not leak out. These sets act as routes for the 
material transport. The approach is based on the Perron-Frobenius transfer opera¬ 
tor, which models how the flow moves distributions of points, instead of individual 
trajectories. The Perron-Frobenius operator is always infinite-dimensional and lin¬ 
ear; the identification of almost-invariant sets is then intimately connected with 
approximating its eigenfunctions. While the Perron-Frobenius operator has been 
a staple of the ergodic and probability theory since the early 20th century, it was 
introduced to the applied, non-probabilistic context by Dellnitz and Junge [11,12] 
as the basis for identification of invariant sets in time-invariant dynamical systems. 
Since then, the theory has been expanded to include detection of almost-invariant 
sets of autonomous systems [18,20], and flows with more general time dependen¬ 
cies [19,21]. 

Spatial invariants of dynamical systems relate to infinite-time averages of func¬ 
tions along Lagrangian trajectories. This connection between the ergodic theory 
and applied mathematics was initially developed in Ref. [43] . Based on these ideas. 
Ref. [56] proposed that even finite-time averages of functions can enable detection 
of geometric structures important for fluid transport, which broadly constitutes the 
mesochronic, i.e., time-averaged, theory of transport in fluids. The utility of time- 
averages has been corroborated on numerical and experimentally-realizable flows 
with simple time dependence [40-42,47,48]. 

The mesochronic techniques have developed in two directions. One focuses on 
computations of ergodic invariant sets using long-time averages of a large set of 
averaged functions. [8, 36, 43, 48] The other, which we follow here, does not aim 
to compute all invariant sets; rather it uses much shorter averages of the velocity 
held itself to identify the character of the deformation, i.e., presence or absence of 
rotation. [46] This is in contrast with mentioned LCS and rotational theories, which 
describe deformation starting from analysis of the magnitude of the deformation. 

Before we dive into calculations, we give a short explanation of the approach. 
Introductory courses in dynamical systems often discuss the stability of a stagnation 
point p in a planar system x = f(x) by looking at its linearization ^ = V/jp • ^ 
around a fixed point p whose stability depends on positions of two eigenvalues of 
the Jacobian V/|p. Instead of computing the eigenvalues, their locations can be 
inferred from the trace and the determinant of V/|p. If the flow is incompressible, 
the trace is zero, so the determinant alone is needed for the full stability analysis. For 
unsteady systems this analysis may not always hold; however. Ref. [46] showed that 
even then similar results can be obtained by looking at the Jacobian matrix of the 
velocity averaged over a Lagrangian trajectory, termed the mesochronic Jacobian 
matrix. Away from fixed points, the calculation does not compute stability, but 
rather the spectral class of the Jacobian: hyperbolic (strain) or elliptic (rotation), 
termed mesochronic classes. Applied to prediction of the oil slick transport in the 
aftermath of the Deep Water Horizon spill [46], mesochronic analysis showed that 
regions of mesohyperbolicity correspond to jets which dispersed the slick, while 
mesoelliptic zones correspond to centers of eddies in which the slick accumulated. 
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The main contribution of this paper is the extension of the mesochronic analysis to 
three-dimensional flows. 

There are several connections of the mesochronic analysis with other approaches. 

• On a fundamental level, averages of functions along trajectories are inti¬ 
mately related to spectral properties of the Koopman operator [35,44,45], 
which is adjoint to the mentioned Perron-Frobenius operator. 

• Greene [26,27] defined the residue criterion in order to predict the order of 
destruction of Kolmogorov-Arnold-Moser (KAM) tori in perturbed Hamil¬ 
tonian maps. The computation of the residue is almost identical to that 
of the mesohyperbolicity indicator for 2D flows. The three-dimensional 
version of the residue criterion [17] also resembles mesochronic indicators 
introduced here. 

• An analysis of oceanic flows based on the Jacobian of the instantaneous, i.e., 
non-averaged, velocity is well-known as Okubo-Weiss-Chong partition [9, 
51,64]; this is the limit of the mesochronic theory as the averaging time 

0 . 

• In the other extreme, as T —>■ oo, real parts of eigenvalues of the mesochronic 
Jacobian relate to Lyapunov exponents and rotation numbers). It can be 
shown under generic conditions, using the polar decomposition of V'i/'T that 
the limit of the eigenvalues of the gradient of the flow map are the Lyapunov 
exponents (see the heuristic discussion in [23], that can be turned into 
rigorous proofs using ergodic-theoretic techniques in [58], see also [1, Chapt. 
3, §9]) 

• Finally, a recent inquiry [15] into connections with Lagrangian Coherent 
Structures, in a form related to an earlier work in Ref. [56], showed that the 
LCS techniques are capable of uncovering some of the boundaries between 
mesochronic classes. 

The paper is organized as follows. In Section 2 we introduce the precise defi¬ 
nitions for dynamical systems we are considering, review the basics of differential 
geometry needed, and re-state the Okubo-Weiss-Chong analysis in these terms. 
Section 3 contains the main result, the 3D mesochronic classification, while Sec¬ 
tion 4 makes connections to Lyapunov, Okubo-Weiss-Chong, and 2D mesochronic 
analyses. In Section 5 we illustrate the technique by a set of analytic and numer¬ 
ical examples; in particular the steady and unsteady Arnold-Beltrami-Childress 
flow [13]. Numerical details are given in the Appendices. The paper closes with the 
discussion in Section 6. 


2. Preliminaries 

Consider a time-varying differential equation 
( 1 ) x = f{t,x) 

with a velocity field / : D C K x —>■ K^. For an initial condition (to,xo) € D 
let 1 1 —> Lp(t^ to, ccq ) denote the solution (or trajectory) of the initial value problem ( 1 ), 
with a;(to) = xq. Throughout the paper we assume for an arbitrary but fixed initial 
time to € duration T > 0 and open set of initial values X{to) C that 
Vcco € A'(to), t !-)• (fi{t,to,xo) exists on the whole interval I := ]to 7 ^o + T], 

For t G I define X(t) := {ip(t,to,xo) G : xq G X{to)} and X := {(t,x) G 
K X : t G /, X G X(t)}. Then X C D hy assumption und for ti G I the map 
(p{tQ -I- T, to, •) : X{to) -G X{ti) is well-defined. In particular, we define the time-T 
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map 

^T-.X{to)^X{to+T), 
iPt{x) ■■= (p{to + T,to,x), 

usually called Poincare map if the equation is periodic, i.e., if for the chosen T, 
f{t,x) = f{t + T,x). 

We are mainly interested in finite-time dynamics for a fixed duration T > 0 
but will also investigate the instantaneous (in the zero-time limit T —>■ O’*") and 
asymptotic (in the infinite-time limit T —>■ -|-oo) dynamics, assuming the solution 
ip{-,to,xo) exists on [fo,oo). 

Observables of (1) are continuous functions F : X ^ K”, which are evaluated 
along arbitrary solutions 1 1 —)■ x{t) on I. They are used to model physical measure¬ 
ments of a state of the system, e.g., the time trace t^F{t, x{t)) might represent the 
ocean temperature recorded by a sensor as it is passively carried by ocean currents 
along the trajectory x(t). A time average or trajectory average Ft : Xfto) —>■ K of 
an observable F on I = [to, to + T], defined by 

^ rto-\-T 

Pt{Xo) ■■= 7^ F{T,x{T))dT, 

is a function that depends on the initial value a;(to) = xq of the trajectory x(t) at 
time to- Trajectory averages depend on the duration T > 0 and can be analyzed 
from the instantaneous, asymptotic, or finite-time perspective. The instantaneous 
case is the most obvious, as lim-r-^o Ft{xo) = F{to,xo), e.g., if F is a component of 
the velocity field itself. Certain choices of F can still provide valuable information, 
as we explain in the next paragraph. Asymptotic analysis studies ergodic averages, 
i.e., limits Ao(a^o) := liniT-^oo A(xo), in case they exist. Ergodic theory analyzes 
limits Foo of observables F which are specified only in general terms, e.g., only 
by the space of functions from which they are drawn. Even in such general cases 
valuable information can be recovered, e.g., on time-invariant measures on the state 
space. [8,48] 

On the other hand, choosing a particular observable can provide us with more 
detailed information. Since the components of the velocity field f(t,x) = [fi{t,x), 
f 2 {t, x), fsit, x)]^ are themselves continuous functions on the time-state space X C 
M X they are observables. One could argue that they are the most distinguished 
observables for analysis of dynamical systems, as they directly provide dynamical 
information about the behavior of the system. We adopt the velocity-as-observable 
viewpoint and analyze the time average of the velocity field, which was termed 
mesochronic velocity [46]. We note that the values of the mesochronic velocity 
appear as quantities of interest in [56] while [46] is the first to look into their 
gradients. 

Definition 2.1 (Mesochronic velocity and mesochronic Jacobian). The mesochro¬ 
nic velocity / : Xit^) C —)■ of (1) on / = [to, to + F] is given by 

frix) := - J f {T,(p{T,to,x))dT. 

The Jacobian matrix V/r containing partial spatial derivatives [V/t]^^ := d[fT]i/dxj 
is termed the mesochronic Jacobian. 
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Note that the spatial derivatives and the averaging over trajectories do not com¬ 
mute, i.e., V/t is not equal to the average of the instantaneous Jacobian over 
trajectories. 

The mesochronic velocity in the instantaneous limit fx - > f coincides with 

the velocity field /. The asymptotic limit for T —>■ oo exists in many cases, for ex¬ 
ample, if the dynamical system is autonomous and volume-preserving on a compact 
domain. 

We use the mesochronic velocity to determine the character of the evolution of 
a material volume (see Section 2.1) by an incompressible dynamical system, which 
satisfies the Liouville condition 

(V • f){t,x) = tiVf{t,x) = 0. 

Although the limits T —)■ 0'*' and T —)■ oo have been studied classically, nei¬ 
ther theory is applicable to transient behavior. On the other hand, the finite¬ 
time analysis of the mesochronic velocity recovers the character of the time-T map 
^ : X{to) X{to + T), which captures transients at the time-scale T [46]. 

2.1. Deformation of a volume cell under a diffeomorphism. We now briefly 
review basic differential geometry that characterizes deformation under a volume¬ 
preserving diffeomorphism -0 using the spectral class of its Jacobian V'0. This 
theory is later applied to time-T maps tpT of dynamical systems over finite time 
intervals. 

Let 'tjj : U ^ V he a diffeomorphism between two open subsets C/ C K^, F C K^, 
with the usual volume measure on K^. We are interested in deformation of an 
infinitesimal volume cell surrounding x £ U as x ^ 4’i^). The central object of our 
interest is the Jacobian matrix Xip ■ U —>■ It is a basic result in differential 

geometry that volumes of a set S C U and its image il’iS) C V are equal if and 
only if |det V'i/’(a;)| = 1- We now restrict our attention to orientation- and volume¬ 
preserving maps Ip, i.e., maps for which det V^/> = 1. 

At the coarsest level, we distinguish between a hyperbolic deformation, when the 
volume cell is deformed along all three spatial dimensions, and the opposite, non- 
hyperbolic character. Let pi, fj, 2 , fJ -3 S C denote the eigenvalues of V'0, assuming 
I Ml I < I M 2 1 < I M 3 1- Different fields of mathematics may interpret presence or absence 
of hyperbolicity differently, e.g., as the material deformation in continuum mechan¬ 
ics, or stability of the map in dynamical systems and control. Since our analysis 
could include both domains, we refer to presence/absence of hyperbolicity, and their 
sub-classification (see below) as the spectral character of the diffeomorphism. 

Definition 2.2 (Hyperbolicity). The map tp is hyperbolic at x if no eigenvalues 
Mi,M 2)M3 of fho Jacobian Xtp{x) lie on the unit circle in the complex plane, i.e., 
V i = 1,2,3, |^i| ^ 1. Otherwise, it is non-hyperbolic at x. In particular, if all 
eigenvalues lie on the unit circle, i.e., V i = 1, 2,3, \pi\ = 1, the non-hyperbolic map 
is elliptic. 

Depending on the complexity of eigenvalues, we further distinguish sellar (saddle- 
like) and helical (spiral/helix-like) character of the deformation under ip.^ 


^From Latin, sella, saddle, and helix, spiral. 
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Definition 2.3 (Sellar and helical deformation). The map tjj is sellar at x if all 
eigenvalues of the Jacobian V'!/'(x) are real and non-defective (their alge¬ 

braic and geometric multiplicities match). If, instead, there is a pair of complex- 
conjugate eigenvalues, the map is helical (spiral-like) at x. 

Deformation at a sellar point exhibits three distinct spatial axes, meeting at x, 
whose directions are preserved under ip. The directions of preserved spatial axes 
correspond to real-valued eigenvectors of Vi/^, while the associated (real) eigenvalues 
/i determine whether the points along the axes are moving away from x, (|/i| > 1), 
moving towards x, (|^| < 1), or remain neutral (|/i| = 1). 

Deformation at a helical point exhibits only a single preserved spatial axis, around 
which a volume cell is rotated, resulting in a single real-valued eigenvalue. As 
Vil){x) is a real matrix, any complex eigenvalues must arise in conjugate pairs, /r, p,. 
The modulus of the complex pair again determines expansion or contraction of the 
material in the rotation plane, while the real and imaginary components of the 
associated eigenvector pair span the rotation plane. 

Since we restrict our analysis to volume-preserving maps ip, existence of an eigen¬ 
value inside the unit circle (contraction) necessarily means that at least one other 
eigenvalue lies outside of the unit circle (expansion), as |detV'0| = |/ri/i 2 /r 3 | = 1. 
All possible combinations are enumerated in Definition 2.4 in which we use the 
symbols -|- and — to denote expansion and contraction directions, respectively. 

Definition 2.4 (Spectral classes). Let ip : U —>■ V he a volume- and orientation¬ 
preserving diffeomorphism, x G U a point, and Vip[x) the Jacobian of ip at x with 
eigenvalues ordered as |/ii| < |/i 2 | < The class of the point x is determined 
according to Table 1 by the number of eigenvalues of Vip{x) that are inside and on 
the unit circle. 

The first four classes are hyperbolic, whereas the remaining cases are non-hyperbolic. 

Informally, we will refer to signatures [-1—[-] and [-[-] of hyperbolic points 

as, respectively, flattening and elongating, due to the shape of a volume cell after 
application of ip, as sketched in Figure 1. 

2.2. Instantaneous deformation by a dynamical system. We now apply the 
classification outlined in Section 2.1 to the time-T map ipr defined in (2) for a 
fixed time interval T to recover the Okubo-Weiss-Chong [9, 51, 64] criterion for 
classification of velocity fields. Instead of forming the classification of ipx based 
on properties of VipT, we chose instead to base our approach on properties of the 
mesochronic Jacobian V/t, resulting in criteria whose expressions are consistent 
across all time scales T. 

Instantaneous classification refers to the infinitesimally small time interval length 
T for which the class of x G under ipT can be inferred from the spectrum of the 
Jacobian Wf{to,x), through the connections with Wipr- Both Jacobians are 3x3 
matrices; for a general such matrix M, the characteristic polynomial is given by 

Pm{s) = det (sId — M) = tr M -I- s tr Cof M + det M 

where 

trM = ^Si, detM = ]^Si, 

i i 

trCof M = = ^(trM^ - {ti Mf) 

Si 2 


( 3 ) 
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Class 

Condition 

[-h +] sellar 

(hyp. saddle) 

Ml, 2,3 G K and 

I/ill < 1 < |/i2| < |/i3|, 

[-h] sellar 

(hyp. saddle) 

Ml, 2,3 G K and 

ImiI < IM 2 I < 1 < IM3|, 

[—1—h] helical 
(hyp. spiral) 

Ml G M, M 2 = M3 ^ and 
ImiI < 1 < lM2i = |M3|, 

[-h] helical 

(hyp. spiral) 

Ml = M 2 ^ K, M3 G ffi and 

ImiI = IM 2 I < 1 < IM3|, 

Neutral saddle 

Ml 2 3^^ ctlld 

ImiI < 1 = IM 2 I < IM 3 I, 

Neutral helix 

Ml = M 2 ^ M3 G K and 

Im*I = 1, 

Pure shear 

Ml = M 2 = ±1, M3 = 1 and 
Vi/) is defective, 


Pure reflection G { 1,1} and 

V'0 is not defective. 

Table 1. Classification of a 3D volume-preserving diffeomor- 
phism, depending on locations of eigenvalues fJ. 1 . 2,3 of the Jaco¬ 
bian matrix (see Definition 2.4). A matrix is defective, or non- 
diagonalizable, when it has less than 3 linearly independent eigen¬ 
vectors. 


The coefficient trCofM is the cofactor trace, also called “second trace”, in [17]; 
intuitively, it is the “variance” of the vector containing eigenvalues. Additionally, 
it follows from (3) that when detM = 1, the second trace is equal to the trace of 
matrix inverse trCofM = = trM“^. Further spectral properties of such 

matrices are summarized in the Appendix A. 

As mentioned in Definition 2.4, spectral class of the time-T map is determined 
by locations of eigenvalues of VipT which, in turn, are roots of the characteristic 
polynomial = det (/rid — VipTix)). 

Expanding Viplx) into a Taylor series around T = 0 yields 

,,, ViPTix)=ld + TVf{to,x) + OiT^) 

(4) 

«Id + TV/(fo,x), 


for T « 0. Consequently 

P,p{n) = det [{fi - l)Id - T\7f{to, x)] 


= T^P, 


- 1 
T 


for T « 0, where P/(A) = det (Aid — V/(to,a^)) = A^ — t/A^ -I- m/A — df is the 
characteristic polynomial of the Jacobian matrix V/ of the instantaneous velocity 
held. In all cases, Jacobian matrices are evaluated at the time-space point {to,x) 
which is suppressed in notation. 
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(A) [- 
mation 


-h] sellar defor- 

(elongating saddle) 



(d) [-h] helical defor¬ 

mation (elongating spiral) 



(b) [-h +] sellar defor¬ 

mation (flattening saddle) 



(e) [-h +] helical defor¬ 

mation (flattening spiral) 



(c) Neutral sellar (neutral 
saddle) 



(f) Neutral helical (neutral 
spiral) 


Figure 1 . Sketches of deformation depending on the spectral class 
of a volume-preserving diffeomorphism i/i at the point x € U (see 
Definition 2.4). Note that the initial and final axes do not need 
to be parallel in general. Pure shear and reflection (including the 
identity) are not sketched. 


Incompressibility of the flow implies 
(5) ^/ = 0 and = 1, 

and therefore the spectral character depends only on the determinant df and the 
sum of minors m/ of the velocity field Jacobian V/. We now state the Okubo- 
Weiss-Chong criteria, using terminology specified in Definition 2.4. 

Theorem 2.5 (Okubo-Weiss-Chong [9]). From the Jacobian of the velocity field 
(velocity gradient tensor) V/(to,a^)> compute its determinant df = det V/(to, a;), 
and cofactor trace mj = trCof V/(fo 5 a:). For T —>■ 0, the point {tQ,x) G I xM.^ is 
hyperbolic if and only if 

df 7^ 0. 

Hyperbolic points can be further classified into four subclasses, depending on signs 
of the determinant df and the quantity 27c?^ -|-4 toj, as listed in Table 2. 

The analogous result for planar dynamics is known as the Okubo-Weiss crite¬ 
rion [51,64]. 


3. Mesochronic classification 

The Okubo-Weiss-Chong criterion (Theorem 2.5) classifies the points in the 
domain of the time-T map ipT according to the spectral character of the Jacobian 
V/ of the velocity field. In general, the spectral character of Vi/’t for T away from 0 
will not be approximated well by the spectral character of V/, except in extremely 
slowly varying flows. In order to capture both deformation at small and large T, 
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df 

27d^j: + 4m^ owe class 

- 

— [-h +] saddle 

- 

-1- [-h +] helix 

+ 

— [-h] saddle 


+ + [-h] helix 


Table 2. Okubo-Weiss-Chong classification based on signs of 
quantities in the first two columns (see Theorem 2.5). 


we replace V/ by the Jacobian V/r of the mesochronic velocity fx, defined as the 
Lagrangian average of the velocity / over the duration T: 

1 

(6) frix) := - J f{T,(p{T,to,x))dT. 

To see why Vfxix) plays an important role in our analysis, we write the solution 
ip{t,tQ,x) in its integral form 

(7) ip(t,to,x) = X+ f{T,ip{T,to,x))dT, 

Jto 

which is a consequence of the fundamental theorem of calculus. The same integral 
appears both in (6) and (7), which we use to write the time-T map as 

(8) ij^rix) = (p{to+ T,to,x) = X+ TfT{x), 
and hence 

(9) VV't(x) = Id + TVfrix). 

Comparing with (4), Viprix) « Id+TV/(to, x) which is valid for T « 0, we see that 
for general time intervals [tg, t^+T] the mesochronic velocity Jacobian V/ 7 ’(a;) plays 
the same role as the velocity field Jacobian V/(fo, x) did for intervals of infinitesimal 
length. 

In analogy to Okubo-Weiss-Chong, we classify a; € X(to) C under the action 
TpT using only information obtained from the mesochronic velocity and its Jacobian 
V/r(x). 

Definition 3.1 (Mesochronic classes). Given a fixed time interval +T], the 

point X G X{to) C is mesohyperbolic if it is hyperbolic with respect to the 
diffeomorphism 'ijjx (time-T map), i.e., no eigenvalues of 'V'lpTix) lie on the unit 
circle in the complex plane. Otherwise, x is non-mesohyperbolic. 

Similarly, the mesochronic class of x is the spectral class of ipT at x, as specified 
by Definition 2.4. 

As explained in Section 2.1, the Liouville incompressibility results in only two, 
instead of three, “axes” in which we understand the spectral class of ipT'- 

(1) number of contracting directions, indicated by labels [—1-+] and [-[-]• 

(2) presence of a rigid rotation, indicated by helix vs. saddle split. 

The four mesohyperbolic classes are formed by choosing an option along each of 
the axes above, with the non-hyperbolic classes acting as boundary cases between 
them. 
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E 

A 

Mesochronic class 

- 

- 

[-h +] mesosellar 

- 

0 

[-h +] mesosellar or shear 

- 


[-h +] mesohelical 

0 

- 

neutral mesosellar 

0 

0 

pure reflection 

0 

-t- 

neutral mesohelical 

-f 

- 

[-h] mesosellar 


0 

[-h] shear or mesosellar 


-f 

[-h] mesohelical 


Table 3. Mesochronic classification based on signs of E and A 
(see Theorem 3.2). The classes with E 7^ 0 are mesohyperbolic. 
The class for E = ±00 is the neutral (non-mesohyperbolic) saddle. 


For conceptual and computational reasons we will formulate quantities E and A, 
each corresponding to one of the “axes” above, such that the signs of their values 
sort the finite-time dynamics around x into mesochronic classes. The number of 
contracting directions will be detected by E, the presence of rotation by A. 

Theorem 3.2 (Mesochronic classification). Let Vfr^x) be the mesochronic Jaco¬ 
bian matrix for the dynamics at the point x G X{to) C and for the time interval 
[to, to + T] o,nd PfW '■= — dp the characteristic polynomial of the 

matrix X/fxix). Define 

dfT^ 

j] _ I _ 

8 - 2m^-r2 - 3dfT^ 

A := - AdjT^"^ - Udjm^T^^ 

- edf-mjT^ 

(ISdjTOj — mP)T^ -\- 

-h (27fi2--h4mpr®. 

If 8 — 2m— SdjT^ = 0, then x is non-mesohyperbolic of the neutral-saddle type. 

If E is finite, the point x is classified into mesochronic classes according to Ta¬ 
ble 3. Mesohyperbolic classes are those for which E is finite and E 7^ 0, i.e., for 
which both 

(11) df ^0 andS- 2m- 3dfT^ ^ 0. 

The distinction between the shear and the boundary saddle class when A = 0 
depends on eigenvectors of the mesochronic Jacobian and cannot be made purely 
based on the spectrum. 

The proof of the theorem will be the result of two lemmas: 
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• Lemma 3.3 defines E and A using coefficients of the characteristic polyno¬ 
mial of the time-T map Jacobian matrix Vi/’t (a;)- 

• Lemma 3.4 establishes relations between characteristic coefficients ofWtpTix) 
and characteristic coefficients of Wfrix). 

Lemma 3.3. Let iPt{x) be the time-T map at a point x G A'(to) C over the 
time interval [to, to + T] and := — t^p? m^p — d^, the characteristic 

polynomial of the Jacobian matrix \'^fT{x). Define 
^ t'tp rn^ 

(12) ty, -I- 771^ -I- 2 ’ 

A = 4(to^ -I- t^) - m^t^ - ISm^t^ -I- 27. 

U tif, 2 = 0, then x is non-mesohyperbolic of neutral-saddle type. 

If E is finite, then the point x is classified into mesochronic classes defined by 
Table 3. Mesohyperbolic classes are those for which E is finite and E yt 0, i.e., 

t., 1 , — m.,p 7 ^ 0 and t.^ -I- -|- 2 7 ^ 0. 


Remark 1. The trace t^ and the determinant d.,p of the Jacobian matrix Vi/’t 
are commonly encountered in matrix analysis due to their simple relationships to 
eigenvalues. To interpret the cofactor trace we can expand it as: 


3 

jjk i=l 


1 


= det Vi/'T • tr 




Under an incompressible flow (detV-0T = 1) the cofactor trace is the trace of 
the inverse of the time-T map Jacobian, as discussed after (3): 


Additionally, we can clarify the meaning of E. Rewrite (12) as 


E = 1 - 2, 

-\-1 

Then through simple algebraic manipulation it follows that 


sgnE = sgn ( \ - 1 ) = sgn(t,/, - to,/,). 

\ ru,/, -|- i / 

Using TO,/, = f,/,-i and rewriting the expressions for f,/, using eigenvalues pi of the 
flow map Jacobian, we obtain that 


SgnE = sgn [{pi - p.^^) -\- {p2 - M2 “ (MiM2 - Mi V2 ^)] ■ 

Proof of Lemma 3.3. The mesochronic class of x is determined by eigenvalues of the 

Jacobian 'S/ipTix), whose characteristic polynomial is P.,j,{p) = m^ — 

with determinant d.ip, trace f,/,, and cofactor trace m.,p of V'0r (see Appendix A). 

The mesochronic class is determined by relative locations of zeros /^^(m*) = 0 
in reference to the unit circle and to the horizontal axis. As mentioned, converting 
the reference unit circle to a vertical axes simplifies the criteria, as they can then be 
computed solely from signs of eigenvalues. To this end, we introduce the conformal 
map r(s) = which maps the left half-plane in C to the inside of the unit circle 
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in C, while preserving the upper half-plane. It follows that the location of zeros 
of the composite function oT with respect to axes is the same as the location 
of zeros of with respect to horizontal axes and the unit circle. Note that the 
inverse r“^(/r) = has a pole at fi* = —1, so no finite zeros in the s-plane will 
correspond to the zero of P^ at n* = —1. For this reason, we will separately treat 
the case when WtpT has an eigenvalue at —1. 

Assuming P^{—1) ^ 0, the composite function is 


with coefficients 


[P^oT]is) = 


nss^ + n2S^ + nis -|- no 


(s-lf 


n3 — 1 TTi-xp d'ip 

Ti2 — 3 iy'ip ~\~ Tn^ “t“ 
ni = —3 -\- t'lj, -\- — 3(i^ 

no = — 1 -f — TTi-^ 


t'(p TTltp 2 

tip TTlip^ 


where the second equalities are obtained by incompressibility dip = Jlti Pi = 1- 
A point X is non-hyperbolic whenever one of the roots of the characteristic poly¬ 
nomial Pip{pt) of is on the unit circle. If Pp,{—1) p 0, this condition implies 

that a purely imaginary number m, a G M is a zero of P^ o F. Substituting into 
the numerator of o F, we obtain the condition nin2 — ngna p 0 for hyperbol- 
icity. This condition translates to a condition on the spectral coefficients of \/tpT- 
tip — rriip p 0. The other case for non-hyperbolicity is Pip{—1) p 0, which by direct 
substitution into Pp, translates to t^ + -I- 2 7^ 0. 

When P^{—1) p 0, we can proceed to a finer classification of the location of 
the roots of Pip. Due to incompressibility, there either have to be two zeros of the 
polynomial Pip inside a circle and one outside, or vice-versa. Under a conformal 
transformation F this condition translates into two zeros of Pip o F having matching 
signs, while the third is of the opposite sign. It follows that the product of zeros s* 
of Pip o F is positive when two zeros are negative, corresponding to two directions of 
contraction under the flow map, while it is negative when there is a single contracting 
and two expanding directions. By using the zeros s* of Pip o F to factorize its 
numerator = ^3ni=i('® ~ ^t) equating the zeroth-order coefficients 

one obtains —ns Jlti = ^o- We therefore define an indicator S as 

^ .t. * * ^0 tip niip 

7 1 * — 51 So So —-— -• 

R-3 tip + TUip + 2 

By this argument, S 7^ 0 implies hyperbolicity: E > 0 implies two directions of 
contraction, one of expansion, while E < 0 implies one direction of contraction, two 
of expansion. When E = 0 one of the roots of Pip lies on the unit circle, which we 
term “neutral”, while the other two directions are contracting and expanding. 

The presence of rotation is indicated by the complexity of zeros of Pip o F, which 
is determined by the discriminant of the numerator, 

(13) VsD = -64(4(m^ - t%) - -I- 27). 

When T> 3 d > 0, all zeros are real and distinct, when T> 3 d < 0, two zeros are 
complex-conjugates of each other, when = 0, one real zero is repeated, while 
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the third is distinct. We therefore define our second indicator A to be 


A := -I?3£,/64 = 4(to^ + t^) - + 27. 

It follows that A < 0 indicates that all eigenvalues of Vi/’t (a^) are real, implying 
that there is no rigid rotation; A > 0 indicates that two eigenvalues are complex, 
so rigid rotation is present. When A = 0 two eigenvalues coincide somewhere along 
the real line. 

Finally, we have to deal with the case when P^(—1) = 0 which is not covered by 
the s-complex plane as defined above, as it implies that the denominator of E is zero, 
i.e., + 2 = 0. If one zero is then the two others have to be /r* = z, 

^3 = —\/z, for some z S C to satisfy incompressibility. Therefore the studied point 
X is non-hyperbolic, with signature [— 0 +]. Notice that if z has an imaginary 
component, its conjugate is z = —1/z, which cannot be, as zz > 0, therefore one 
zero at —1 implies that the other two zeros are necessarily real. It follows that 
P,0(—1) = 0 indicates a neutral saddle case, even if E cannot be evaluated. □ 


We now relate the characteristic polynomials of VV't and V/t. 

Lemma 3.4. The characteristic polynomial of the Jacobian matrix VV't O'nd the 
characteristic polynomial of the mesochronic Jacobian V/t are given by the expres¬ 
sions 

M t^fX + TTl'ipp d.0, 

(14) , , , 9 

Pj(A) = — tjfX^ mjX — dj. 

The coefficients are linked by the expressions 

t^ = 3 + Tt j 

(15) = 3 + 2Tt^ + 

= 1 + Tt^ + + T^dp 

where T is the length of the averaging period in /t. 

Moreover, the incompressibility condition = 1 imposes the relation 

(16) + Ttoj + T^d^ = 0 

on the coefficients of for T 


Proof. The connection between the two Jacobian matrices V'0t and V/t is given 
by relation (9): V^’t = Id + TV/t. The characteristic polynomial of V^/t can 
be re-written using the characteristic polynomial P^ of V/t 




det(^Id — V^t) = det 


(^^W-V/t) 



Using the notation for coefficients of Pj from (14), we can expand the last line to 
obtain 

P^(/r) = - {3Ttj)fi^ 

+ (3 + 2Tt^ + T'^mj)p 

— (1 + Tty + T'^rrij + T^dj). 
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By comparing coefficients with the general expression for P^, the statement of the 
theorem follows. The incompressibility condition from the statement of the theorem 
is a consequence of substituting = 1 for any T > 0. □ 

We now combine the preceding Lemmas to give the proof of Theorem 3.2. 

Proof of Theorem 3.2. We start with the expression for the P.^{—1) ^ 0 condition, 
which was stated as + 2 ^ 0. Using (15) we can rewrite the condition as 

8 + + mjT^ 7 ^ 0 . 

Since tj, rup and are related through the incompressibility constraint t^+m^T+ 
djT^ = 0 derived in (16), we can formulate the condition in two alternative ways: 

8 + 2tfT - djT^ 7 ^ 0 

8 - 2m 7 ^ 0 . 

The other mesohyperbolicity condition — m,/, 7 ^ 0 translates into either 
8 dyT^ 7 ^ 0 or 8 T(t^ + TOjT) 7 I 0 . 

Reformulating the expressions for E and A in Lemma 3.3 

+ m^p + 2 tj/j + m.p, + 2 

A = 4(m^ + 1^) - + 27, 

using dp mj and through (15) we obtain 

tf+ Tmf 
yj _ _ 'j' _ I _ J 

8 + iTtj + T'^mj 

A = (4m^ — + (18TO^tj — 

+ (27m| + 18 to^t 2- - 4tpT^ 

+ + 27tjT‘^. 

Since tp mj, and d^ are related through the incompressibility constraint (16), 
we can re-formulate the expressions using either or 

dfT^ dfT^ 

V _ _£_ _ _£_ 

“ 8 -h 2tfT - dfT^ “ 8 - 2m^-T2 - ’ 

A = - -h edjtfT'^ 

+ {27df + 2tj)djT^ + 

= - 4djT^'^ - UdjmfT^^ 

— 13djTOjT^° — hd^m^T® -f (18d^mj — 

+ ISdjmjT'^ + (27d®- -h 4mj)T^. 

To emphasize the connection to the OWC expressions (Theorem 2.5), we will use 
the dy and versions of the formulas. 

The statement of the proof is equivalent to Lemma 3.3 where the criteria are 
expressed in terms of the spectral coefficients of V/r instead of V' 0 t- D 
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In summary, to identify the mesochronic class of a point x, take the following 
steps: 

(1) compute the Jacobian Vfrix) (details in Appendices B and C), 

(2) evaluate A and E using (10), 

(3) use Table 3 to identify the mesochronic class. 

4. Limits and boundary cases 

When, respectively, T —)■ 0 and T —>■ oo, we show that mesochronic classifi¬ 
cation limits to classical Okubo-Weiss-Chong and Lyapunov exponent analyses. 
Additionally, if the dynamics evolves on an invariant plane in 3D space, certain 3D 
mesochronic classes have their counterparts in the 2D mesochronic classification [46] . 

4.1. Instantaneous limit: The Okubo Weiss Chong criterion. The rela¬ 
tion between flow map and mesochronic velocity (8) can be rewritten as fxix) = 
ip(to+T^o,x)-x ^ ^ consequence, the averaged field /t tends to / pointwise as 

T ^ 0'^ and V/ —V/. By continuity of eigenvalues with respect to matrix 
elements, >■ d/, xn/, and t^^tf as T^0~^. As all three of these quan¬ 

tities are finite, the mesochronic incompressibility criterion (16) is trivially satisfied 
in the limit. Additionally, due to incompressibility of the vector field, it holds that 


Theorem 4.1. Suppose that at the point (to,x) the differential equation (1) is in¬ 
stantaneously hyperbolic by the Okubo-Weiss-Chong (OWC) criterion. Then, there 
exists fynin > 0 such that the point {to,x) is also mesohyperbolic with respect to all 
time intervals [fy, fy + T] for which T < fynm. 

Proof. According to Theorem 2.5, we obtain that df fy 0. Since the maps T i—>■ V/t, 
T I—> dy^, and T are continuous, the instantaneous hyperbolicity will imply 

mesohyperbolicity for some small Tmin > 0. The continuity means that there exists 
dmin > 0 on which continuity of the maps T i—)■ dj(T) and T i-A- "iff d —8 

implies 

for T G [0,Tmin]. By virtue of Theorem 3.2, the point (to,x) is also mesohyperbolic 
with respect to [0, T] and the proof is complete. □ 

As a consequence, the signs of the indicators A and E (10) have the following 
limits 

^ J 

sgn Zj ->■ sgn a /, 

sgnA sgn(27d^ -\-4mf). 

By comparing mesochronic classification criteria (Table 3) in the limit T ^ 0 with 
the instantaneous OWC criterion (Theorem 2.5), we can conclude that mesochronic 
classification reduces to OWC classification in the T —>■ 0 limit. Put differently, 
mesochronic classes generalize OWC classes to time intervals [fy, fy + T] with T > 0. 

4.2. Asymptotic limits. For autonomous dynamical systems defined over com¬ 
pact domains, asymptotic rates of deformation are defined by Lyapunov exponents. 
The finite-time analogs, termed Finite-Time Lyapunov Exponents (FTLE) [28,60] 
are defined using the polar decomposition of the Jacobian matrix of the time- 
T map V'T- Let V'0 t = be the polar decomposition of Vi/'t where 
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= WiI’tWi/jt- Eigenvalues of [Vi/'tI are non-negative, singular values of 
Vi/'T) so we can define Finite-Time Lyapunov Exponents ai to be the exponential 
growth/decay rates of the singular values, i.e., we represent the singular values of 
X/i/jt as 

It is a well-known fact in matrix analysis that when a matrix is normal, i.e., 
unitarily diagonalizable, its singular values are equal to absolute values of its eigen¬ 
values. In the language of this paper, this means that when Vi/’t is normal, positions 
of the eigenvalues of Vi/'t in reference to the unit circle are determined by the 
signs of the Finite-Time Lyapunov Exponents ai- It follows that E > 0 (Table 3) 
implies that both two eigenvalues of V' 0 t are outside of the unit circle and that 
two Finite-Time Lyapunov Exponents are positive, when is normal. 

Unfortunately, V'0r is not generally normal for any finite T, i.e., its eigenvectors 
are not orthonormal. However, Ref. [23] shows, although non-rigorously, that for 
smooth ergodic systems which have distinct Lyapunov exponents both real parts 
of eigenvalues and singular values of W'tpT can be written as as T —)■ oo. As 
the sign of E indicates the number of eigenvalues of Vi/’t outside the unit circle, 
which, in turn, is determined by real parts of logarithms of those eigenvalues, we 
conclude that, when T —)■ oo, the sign of E indicates whether two or one Lyapunov 
exponents are positive, assuming that the conjecture in Ref. [23] holds. 

4.3. Recovering the 2D mesochronic deformation criterion. The supple¬ 
ment to Ref. [46] presented a derivation of the criteria for mesohyperbolicity for 
2D (planar) differential equations. Planar differential equations can be trivially 
embedded into a 3D state space by adding a third state with trivial (zero) dynam¬ 
ics. We use such an embedding to demonstrate how the 3D mesochronic criteria 
(Theorem 3.2) specialize to the 2D criterion. 

Let g : / X —>■ be a incompressible (V • g = 0) vector field, with 

mesochronic Jacobian ’S/gx^x) for x G K^, with the spectrum ag = {Ai,A 2 }, and 
trace and determinant tg = Xi + X 2 , dg = A 1 A 2 . The incompressibility implies 

(17) tg + Tdg = 0, 
with eigenvalues then given by 

( 18 ) Xx2 = -^dg±yiT^dg-i)d~g. 

Ref. [46] studied only dg ^ 0, noting that dg = 0 results in Ai = A 2 = 0. The 
time-T map of the velocity field x = g{x) is hyperbolic at Xq if it preserves two 
distinct real spatial axes, which is analogous to the definition in Section 2.1. Con¬ 
sequently, Definition 3.1 retains its meaning for the 2D case. The discriminant 
7 X 20 ■= {T^dg — A)dg (cf. 3D discriminant (13)) indicates mesohyperbolicity if 
7 X 20 > 0 and mesoellipticity otherwise. 

We embed the vector field g in the 3D state space by defining a 3D vector field 
for 7 ; G as 

:= 

V/»= 

where we take Zi _2 •= [^ 1 , 22 ]^ to be the first two components of the vector z = 
[zi, Z 2 , 23 ]^. Eigenvalues of the Jacobian of the 3D averaged field are = {Ai, A 2 , 
A 3 }, with A 3 = 0. The spectral quantities m^, and dj then reduce to analogous 
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quantities of 


tf — 

= Ei=l = dg, 

= 0 . 

As a consequence, the 3D incompressibility condition (16) reduces to the 2D incom¬ 
pressibility condition (17). 

The 3D conditions for mesohyperbolicity (11) reduce to 
dj 7 ^ 0 and 4 — dgT"^ ^ 0. 

Since = 0 as noted above, it follows that the flow is not 3D-mesohyperbolic, as 
it is expected, as the third coordinate is always preserved due to the construction 
of the 3D flow. 

The indicators E and A evaluate to 

E = 0, A = -mjT^irrifT - 4) = -d?r®X> 2 D- 

Therefore, the sign of the A indicator is determined by the 2D mesohyperbolicity 
criterion. If A > 0, T> 2 d < 0 and the two eigenvalues are non-real and lie on the unit 
circle due to incompressibility constraints. If A < 0, I? 2 D > 0 and the eigenvalues 
are real, one is contracting and the other expanding. 

For planar dynamics, I? 2 D = 0 implies that the dynamics are either a pure shear, 
when the geometric multiplicity of the eigenvalue is 1 , or trivial, i.e., tprix) = x. 

Incompressibility in the planar case is more restrictive than in 3D, yielding only 
two structurally stable cases: mesohyperbolic Ai G M and mesoelliptic \i G C, 
which intersect at the pure reflection/shear case Ai = A 2 = 1. The derivation of 
the mesohyperbolicity criterion for the 2D case therefore relied only on detection 
of real vs. complex eigenvalues, which is the reason why I? 2 D in (18) is taken as 
the sole 2D mesohyperbolicity criterion, without resorting to a more complicated 
calculation. 


5. Examples 

5.1. Linear time-invariant velocity fields. In Table 4 we compute explicitly the 
values of the indicators E and A for a simple class of linear time-invariant systems 
whose V' 0 T is constant, and given by the polar decomposition: 

/ cos ujT —sincuT \ / 

Ipj' = [ sina;T cosujT I [ 

\ 1 / V 

In this parametrization, we can independently manipulate rates of strain Ai^ 2,3 as 
well as the rate of rotation uj present in the system. As two of the rates have the 
same sign, unless one of the rates is zero, we choose to order the directions by 
setting sgn Ai = sgn A 2 , which means that the third direction is of the opposite sign 
A 3 = — Ai — A 2 , due to incompressibility. While these systems do not represent a 
broad range of dynamical systems, we have a good understanding of their dynamics 
so it is instructive to see how their properties are reflected in the mesochronic 
classification. 

First, when all rates Ai are non-zero, all points are mesohyperbolic as E is con¬ 
stant and non-zero; presence or absence of rotation determines whether a point is 
a (mesohyperbolic) saddle (uj = 0) or a helix (cu 7 ^ 0). The signature [-h] or 
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V’T 


/ cosujT — sinojT \ / 

= I sinujT cosojT 1 exp I A 2 T 
\ 1/ V -(Ai+A2)t 

(a) Form of the considered LTI systems 


U) 

Ai, A 2 

S(T) 

A(T) 

0 

Ai ■ A 2 > 0 

— tanh AiT/2 

X tanh X 2 T /2 

X tanh(Ai + \ 2 )Tl‘l 

— 64sinh^(Ai — A 2 )r/ 2 x 
sinh^(2Ai + A 2 )r/ 2 x 
sinh2(Ai -|-2A2)r/2 

0 

A = Ai = A 2 0 

tanh AT — 2tanhAT/2 

0 

0 

Ai = A 7 ^ A 2 = 0 

0 

-64smh^{Ar/2)x 

sinh 2 (Ar) 

7^0 

A = Ai = A 2 7 ^ 0 

tanh( AT) (cos(a;T) — cosh( AT)) 

16sin^ {ujT)x 

cosh(AT)+cos(tjT) 

[cos(a;T) — cosh(3AT)]^ 

7^0 

Ai = A 2 = 0 

0 

64 sm‘^{uiT/2)sm^{ojT) 



(b) Values of A and E 


U) 

Ai, A 2 

sgnS(r) sgnA(T) 

Mesochronic class 

0 

Ai ■ A 2 > 0 

- sgn Ai - 

[-h +] saddle (Ai > 0) 

[-hi saddle (Ai < 0) 

0 

A = Ai = A 2 7 ^ 0 

— Sgn A 0 

[-h +] saddle (Ai > 0) 

[-hi saddle (Ai < 0 ) 

0 

Ai 7 ^ A 2 = 0 

0 

neutral saddle 

2D mesohyperbolic 

7^0 

A = Ai = A 2 7 ^ 0 

— sgn A - 1 - 

[-h +] helix (Ai > 0) 

[-hi helix (Ai < 0) 

7^0 

Ai = A 2 = 0 

0 -L 

(c) Mesochronic classes 

neutral helix 

2D mesoelliptic 


Table 4. Mesochronic classes for linear time-invariant (LTI) sys¬ 
tems of the form (a). Values of signs in (c) hold generically, except 
on a non-dense set of periods T where they are zero, as determined 


by values uj and A. 


[—h +] of the saddle is then determined by the sign of the pair Ai, 2 - When the two 
rates match exactly, Ai = A 2 , it implies A = 0. Since the quantities A and S are 
functions of the spectrum, they alone are not enough to detect whether associated 
directions align (shear) or not (saddle). In the case of the systems derived, we know 
that those two directions correspond to independent eigenvectors, which means that 
the point is a saddle. 

If one of the rates is equal to 0, it always implies that S = 0, which is classified 
as one of the neutral 3D mesochronic classes. In that case, A corresponds to the 2D 
mesohyperbolicity indicator T> 2 d, as described in Section 4.3. Again, the presence 
a; yf 0 or the absence of rotation w = 0 is reflected on the sign of A, where A > 0 
corresponds to the former, and A < 0 to the latter case. 
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The magnitudes of A and S grow exponentially in most of the cases; however, no¬ 
tice that in the presence of rotation, a periodic function multiplies the exponentially- 
growing magnitude, resulting in A(T) = 0 periodically. This means that there is a 
potential for resonance, i.e., if T is a multiple of the period of oscillation, dynamics 

momentarily appears to be on the boundary behavior between [—1-+] and [-h] 

saddle mesohyperbolicity, or even a pure reflection when A = S = 0. This choice 
of T is, of course, highly unlikely without a prior knowledge of w. 

In summary, analysis of simple linear systems shows that mesochronic classes 
correctly reflect our intuition about presence of stretching and rotation in linear, 
time-invariant flows. 


5.2. Arnold—Beltrami—Childress Flow. The Arnold-Beltrami-Childress (ABC) 
flow [13] is a kinematic model of an incompressible fluid flow evolving in a three- 
dimensional periodic domain. Even though the system of ODEs specifying the ABC 
flow is simple, it exhibits a variety of different behaviors and has been used as a 
test-bed for various computational algorithms [7,8,20,28]. 

The ABC flow evolves on a 3-torus in periodized state variables {x, y, z) G 
[0, 27r]'^ = T^. Dynamics depend on parameters A, B, C, L> G R and are specified by 
differential equations 

i: = A(t) sinz-I- Ccosy 
(19) y= B sin X + A{t) cos z 

z= Csmy+ Bcosx, 


where the time-varying parameter A(t) is given by 

A{t) = A + Dt sint. 


If D = 0, the equations are autonomous; if, additionally, any other parameter is 0, 
the system is integrable. [13] 

The linearization along a solution p{t) = {x{t),y{t), z{t)) of (19) is given by 


( 20 ) 




0 —Csin.i/(t) A{t) cos z(t) 

Bcosx{t) 0 — A(t)sin 2 (t) 

— Bsina;(t) Ccosy{t) 0 




V f(x,y,z) 


The determinant and the sum of minors are given by the expressions 


df = det V/ = A{t)BC{ cos x cos y cos z 

— sin X sin y sin z) 

m/=trCofV/= A(t)Bsinxcosz 
+BC sin y cos X 
-|-A(<)C'cosy sin z. 

These expressions can be used to evaluate the OWC criterion for the instantaneous 
hyperbolicity according to Theorem 2.5. 

5.2.1. Integrable case. We briefly discuss the case A = D = 0 ^ A(t) = 0 analyti¬ 
cally. From (19), we derive that 

z = Cy cos y — Bx sinx = BC sin x cos y — BC sin x cos y = 0. 

Thus, for the initial condition p(0) = (xq, yo, zq), z(t) = zq + {C sinyo -|- B cosxo)t 
for all t G [0, T], 
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Furthermore, if we write a := x + y, S := x — y, then the ABC flow can be 
rewritten as a decoupled second order system 


a = BC cos a 
S = — BC cos S 
z = 0, 

which is a direct product of two pendulum-like equations. Solutions of the first two 
components can be written in terms of integrals of Jacobi elliptic functions, and 
it follows that the system is integrable. A similar argument follows in case when 
either A, B, or C are zero, in addition to D = 0. 


Lemma 5.1. All points {x,y,z) in the state spaee of the system (19) with Aft) = 0 
are non-mesohyperbolic over any interval + T\. 


Proof. When A{t) = 0, the matrix defining the linear system of equations (20) 


is block-diagonal, with blocks 


and 0 on the diagonal. The 


0 —C sin y{t) 

B cos x(t) 0 

fundamental matrix of the system is, therefore, also block diagonal, with value 1 
on the diagonal corresponding to the exponential of the block 0 in (20). Since 1 is 
then in the spectrum of the time-T map Jacobian, {x,y,z) is non-mesohyperbolic 
for all T. □ 


Remark 2. If Aft) = 0, the mesochronic class of a point (x, y, z) does not depend on 
the value of z since the Jacobian matrix in (20) does not depend on the z-coordinate 
of the solution around which we linearized. 


5.2.2. Steady non-integrable case. The structure of the invariant sets in the state 
space of the ABC flow for parameters A = -s/j, B = C = 1, ZJ = 0 is well 
studied analytically [13] and numerically [8]. The state space contains six interwoven 
vortices with the space between them filled by chaotic dynamics (Figure 2). We 
place a grid of 400 x 400 initial conditions on the (x, y) face of the periodicity cell 
and calculate tj, mj, dj for time intervals of different lengths. Other details about 
numerics are given in Appendix C. 



(a) Isometric view (periodicity cube (b) Slice through z = 0. 

rescaled to unit sides). 


Figure 2. Invariant sets in the state space of the Arnold- 
Beltrami-Childress flow at A = i/S, B — v^, (7 = 1. Regular 
vortices are colored, the space between them is the chaotic zone. 


To give a sense of time scales involved in the system. Figure 3 shows several 
trajectories (pathlines) within a single vortex, simulated for various durations T. 
Trajectories inside vortices take approximately T = 3 to cross one periodicity cell. 
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The two other panels in Figure 3 show that the vortex rotates around its axis while 
the inner layers move at slightly faster speeds than its outer layers. 



(a) Pathlines for T = 1, 2.5, 5 
in a periodized unit cell 



(b) Advection of points for T = 
0,0.75,1.50,..., 4.50 (color in¬ 
dicates time). 





y X 

(c) Advection of points for T = 0,1.5, 3.0,4.5 (color indicates 
time). 


Figure 3. Pathlines and advection patterns for the steady ABC 
flow (19) with A = -y/S, B = C = 1, D = 0. Initial conditions 
clouds were sampled at t = 0 uniformly from the z = 0 section of 
the two vortices in Fig. 2b; panel (b) samples the vortex centered at 
X = Ttt/S, y = 7r/2; panel (c) samples the vortex with the central 
axis along the line x = -k 12. 

To detect non-mesohyperbolic behavior we need to numerically evaluate condi¬ 
tions (3.3) 

(21) ^-^-1=0^ or -I-2 = 0, 

written here using the traces of the flow map and the relation = ty,-i to estimate 
their growth more clearly. These conditions are difficult to reliably compute in the 
face of numerical errors that will almost-certainly result in non-zero quantities. 

Numerically, we determine when the criteria are satisfied using a numerical tol¬ 
erance determined by estimating the growth rate of and with increasing 
length T of time intervals. We infer the behavior of these quantities in linear time- 
invariant systems in order to establish the correct baseline — the criterion has to 
reflect correctly our knowledge of linear systems if we are to even consider using it 
for classification of nonlinear systems. In what follows, we derive an empirical es¬ 
timate of mesohyperbolicity, using a quantity termed numerical hyperbolicity. Our 
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considerations will further rely on T —>■ oo arguments as our intent is to estimate 
the behavior of the flow beyond the time interval in which we are sampling it. 

For a linear, time-invariant system, eigenvalues of the time-T map are given 
either by or by ±e“''‘^. 

In both cases, as T —>■ oo, 

U + I ^maxi AiT 1 4 1 -/• 1 ol ^max^ AiT 

I ^ ' "I ^ 

To account for exponential growth, we set the numerical tolerance of mesohyper- 
bolicity based on logarithms of expressions (21) 

hi ■■= ^ log \t^ - I , := ^ log |4 -h + 2| . 

(In all expressions we omit dependence on state variables and time interval for 
shortness). 

Non-zero values of either hi or /12 are signs of mesohyperbolic behavior; con¬ 
versely, we need either one of them to be small to declare non-mesohyperbolicity. 
In nonlinear flows, we do not expect that hi ^2 will be entirely independent of the 
value of T. Nevertheless, in ergodic regions, [67] we expect convergence in mean as 
T —)■ c». 

Even the rate of convergence to the mean is not uniform: in regular ergodic 
regions, e.g., vortices of the ABC flow, the expected decay is in strongly 

mixing regions, conjectured to be embedded within the chaotic region, the expected 
decay is i.e., similar to the Central Limit Theorem for i.i.d. random vari¬ 

ables. Initial conditions that are neither regular nor strongly mixing may potentially 
have an even slower decay of variance [55], T~°‘ for any 0 < a < 1/2. Values of 
h at those points would then still grow as ~ volume of weakly 

mixing zones is small in systems containing Kolmogorov-Arnold-Moser-type dy¬ 
namics, [8, 54, 63] and therefore we do not expect those values to occur as major 
features in the histogram of h. 

A good quantitative criterion for deciding whether a point is mesohyperbolic 
should estimate whether the smaller of the two hi ^2 


max |Ai| ~ min{/ii, h 2 } —>■ 0 

i 


is “sufficiently” close to zero. Under a conjecture that some variant of the Central 
Limit Theorem (CLT) holds for the estimate of max^ jA^j we can test whether the de¬ 
viation of our estimator min{ft,i, /i 2 } from the hypothesis of non-mesohyperbolicity 
maxi [Ail = 0 is normally distributed, i.e., whether numerical mesohyperbolicity h, 
defined by 


( 22 ) 


h = |min{/ii, /i 2 }| Vt, where, as before 


hi ■■= y log 


dfT^ 


h 2 ■■= y log 


8 - 2m 


is small, h < e, for some (small) constant e. When this is indeed the case, then we 
are reasonably confident that with longer T the estimated eigenvalues of the flow 
map would indeed have a unit modulus and we empirically declare that the point 
is not mesohyperbolic, as classified by Theorem 3.2. 

Is there a fixed value of e that could help us decide whether the hypothesis of 
mesohyperbolicity for a point holds? To answer this question standard CLT re¬ 
sults from probability theory, e.g., Lindeberg-Feller [14, Thm. 3.4.5], would employ 
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higher order moments of the distribution of samples of the random processes hi 2 - 
Such a constant £ would then turn our empirical criterion into a statistical test, 
where h would be the z-score for testing the hypothesis of whether a point under 
consideration is mesohyperbolic or not. Unfortunately, we cannot assume to know 
how the higher moments behave, despite existing work on CLTs in the context of 
dynamical systems [24,25,57,66,67]. 




logiol'j) 


(a) Steady ABC flow 


(b) Unsteady ABC flow 


Figure 4. Distribution of numerical mesohyperbolicity (22) in 
steady and unsteady Arnold-Beltrami-Childress flows for differ¬ 
ent lengths T of integration intervals, computed on a uniform 
grid of 400 x 400 points. Mesohyperbolicity was declared for 
h > 10“°'^ = 0.5 (dotted vertical line). 


In lieu of rigorous results, we choose to proceed empirically and set the cutoff 
value e based on distributions of the numerical mesohyperbolicity h. Figure 4a 
shows histograms of numerical mesohyperbolicity h for a range of values of T, con¬ 
forming well to expectations. As T increases, the distribution of h changes from 
a fairly flat distribution (T = 1) to a bimodal distribution with well-separated 
peaks. Figure 5 shows that each mode of distribution of h corresponds, respec¬ 
tively, to vortices and to the large chaotic region between them as T —>■ oo. Based 
on these results, we declare numerical non-mesohyperbolicity using the cutoff pa¬ 
rameter value 


h < e, with £ = 10 

The Okubo-Weiss-Chong criterion requires df ^ 0 for non-hyperbolic sets. For 
the given parameters at z = 0, 

df(x, y, z) = ABC {cos x cosy cos z — sin cc sin y sin z) 

= VGcosxcosy. 


Therefore, a solution (/?(•, 0,p) with p = (x,y,0) is instantaneously hyperbolic ev¬ 
erywhere except along the lines 


(23) 


TT StT TT StT 



The mesochronic partition for T « 0 is illustrated in Figure 6, and due to short 
integration period T, matches exactly the Okubo-Weiss-Chong partition. Notice 
that the conventional intuition about vortices being “elliptic” structures cannot be 
inferred from short integration times, as for T « 0 almost the entire space is me¬ 
sohyperbolic (non-mesohyperbolic lines (23) are difficult to sample numerically). 
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(d) T = 50 


Figure 5. Spatial distribution of numerical mesohyperbolicity on 
a plane in the state space of the steady Arnold-Beltrami-Childress 
flow. Color is logj^Q h with h defined as in (22). 
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Figure 6. Mesochronic classes in the x — y plane of 2 ; = 0 slice 
of the Arnold-Beltrami-Childress flow state space for T = 10“^. 
As T « 0, the mesochronic partition is virtually identical to the 
Okubo-Weiss-Chong partition. 


Nevertheless, neutral mesohelical regions roughly coincide with locations of vor¬ 
tices, while the chaotic region between vortices contains a mixture of all four classes 
of mesohyperbolicity. Comparing with Figure 2, we see that the boundaries of 
mesochronic classes do not align with boundaries of invariant sets. 

Increasing T results in the sequence of images shown in Figure 7, where the 
mesohyperbolic regions are shown in the left column, and the non-mesohyperbolic 
regions in the right column. As the averaging period is increased to T = 1, partitions 
deform, but remain largely uncorrelated with invariant features. As we increase 
T beyond 1, non-hyperbolic behavior significantly re-appears along the interface 
between different mesohyperbolic classes. Parts of boundaries of mesohyperbolic 
zones start to align with invariant vortices. Since level sets of any function averaged 
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(a) Mesohyperbolic (b) Non-mesohyperbolic 


(c) FTLE 


Figure 7. Distribution of mesochronic classes and the Finite- 
Time Lyapunov Exponents (FTLE) field on the x — y plane oi z = 0 
slice of the steady ABC flow for times (in rows) T = 1, 5, 10, 50. 


for a sufficiently long time will partition the state space into invariant sets [8,48], this 
is expected. Notice that the non-mesohyperbolic regions appear almost exclusively 
inside invariant vortices. As T is increased even further, the non-mesohyperbolic 
zones grow inside the invariant vortices and eventually completely occupy them. In 
the chaotic zone, we see disappearance of mesohelical dynamics, with only saddle 
mesohyperbolicity remaining, which matches the asymptotic analysis in Section 4.2. 
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5.2.3. Unsteady case. We now keep the parameters A = ^/3, B = C = 1 as 
before, but set D = \, which results in the unsteady variation in the coefficient 
A{t) = + tsint. The unsteady ABC flow has not received as much analytic 

attention as the steady case; however, it was used as a demonstration of numerical 
techniques for computation of the flow map Jacobian Vi/’t in [7]. Figure 8 shows the 
same sets of initial conditions used in demonstrating the steady flow (Figure 3), but 
now advected by the unsteady flow for comparison. Notice that the initial advection 
patterns are similar, until A(t) starts significantly deviating from its constant term 
A. 



(a) Pathlines for T = 1, 2.5, 5 (b) Advection of points for T = 

in a periodized unit cell 0, 0.75, 1.50,..., 4.50 (color in¬ 

dicates time). 




(c) Advection of points for T = 0, 1.5, 3.0, 4.5 (color indicates 
time). 


Figure 8. Pathlines and advection patterns for the unsteady ABC 
flow (19) with A = -\/3, B = v^, C = I, D — 1. Initial conditions 
clouds were sampled at f = 0 uniformly from the z = 0 section of 
the two vortices in Fig. 2b; panel (b) samples the vortex centered at 
X = Ttt/S, y = 7r/2; panel (c) samples the vortex with the central 
axis along the line x = 7r/2. 

Figure 4b shows the histogram of numerical mesohyperbolicity while Figure 9 
shows the spatial distributions of (non-)mesohyperbolicity classes, determined using 
the same numerical criterion as in the steady case (22). For short intervals T = 
0.5, 1.0, mesochronic classification of the flow is similar to the steady case. This 
is expected as the magnitude of A{t) is dominated by the steady component. As 
time evolves, non-mesohyperbolic regions in the flow are destroyed, and the obvious 
split between two behaviors that was observed in the steady flow (Figure 7) is not 
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present here. Remnants of the axial vortex in the left and two “eyes” of vortices in 
the right sides of images are visible both in mesochronic and FTLE partitions. 
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(c) FTLE 


Figure 9. Distribution of mesochronic classes and the Finite- 
Time Lyapunov Exponents (FTLE) field on the x — y plane of z = 0 
slice of the unsteady ABC flow for times (in rows) T = 1, 5, 10, 50. 

Figure 10 illustrates transport of initial conditions sampled from several regions 
in the state space of the unsteady ABC flow. The first two rows show advection 
from patches chosen as subsets of regions that are mesohyperbolic for integration 
times T = 0.5,1, 5. The last row shows advection from a patch that straddles several 
mesochronic regions at T = 5. 

Advection up to t = 5 demonstrates that the initial conditions selected from 
single mesochronic regions (central column, top two frames) do not disperse much, 
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(a) Initial patch (hatched) at t = to = 0 (b) Tracer at t = 5 


(c) Tracer at t = 7 


Figure 10. Material advection in the unsteady ABC flow. Rows 
show clouds of 10^ points from regions that are, respectively, meso- 

helical [-h], [-h] and mixed-mesochronic-class between t = 0 

and t = 5. The first column shows the selection patch at time t = 0, 
overlaid from the lower-left square of the third (T = 5) panel in 
Fig. 9a, with two additional times shown in the adjoining columns. 
Point clouds are graphed as they project onto z = 0 plane. 


i.e., stay coherent. Initial conditions from a mixture of regions show considerably 
more dispersion. 

Advection up to t = 7, shown in the third column, demonstrates how the patches 
of initial conditions evolve beyond the interval T = 5 that was used to generate 
mesochronic classes. All material patches at this point show considerable growth; 
arguably, the patch in the last row again shows the largest dispersion. 

We conclude this section with Figure 11 showing images of orbits, i.e., pathlines, 
of material points advected by the unsteady ABC flow through space over 5 time 
units. Two mesohelical patches used to initialize points are the same as in Figure 10; 
the two additional mesosellar patches were initialized similarly. The two mesosellar 
sets are sampled from a narrower region, so it is not surprising that they remain 
more tightly packed than mesohelical sets. We believe that the major distinction 
between the top and the bottom row is in the internal order of orbits within each 
bundle. Mesohelical sets seem to preserve the internal order, i.e., despite the torsion 
of the bulk, one can still observe an ordered rainbow of colors toward the end of the 
orbits. On the other hand, both mesosellar sets show more internal mixing of colors. 
While this is potentially a circumstantial observation here, it may be an interesting 
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(a) [-h] mesohelical patch 


(b) [-h +] mesohelical patch 




(c) [-h] mesosellar patch 



(d) [-h +] mesosellar patch 



Figure 11. Orbits (pathlines) of 100 points initialized from two 
mesohelical patches shown in Figure 10 and two additional mesosel¬ 
lar patches, advected for time t G [0, 5]. Color is added to illustrate 
internal ordering of the material. 


point to explore in the future, due to the association of hyperbolic saddles with 
mixing and bulk rotation with lack of mixing. 

6. Discussion 

Mesochronic analysis is a local analysis in that it classifies a single point based 
on the properties of the local deformation gradients along trajectory emanating 
from it. Nevertheless, the hope is that sets of initial conditions selected using 
the local mesochronic criteria will collectively stay coherent on a macroscopic level 
and, potentially, deform as a bulk in the way suggested by the local analysis. The 
contributions of this paper are in extension of two-dimensional theory of [46] to 
three-dimensional flows. 

We have extended the concepts of mesohyperbolicity (local strain over finite 
times) and mesoellipticity (local rotation over finite times) to three dimensions, 
which allows for co-existence of the rotation and the strain. In turn, we differentiate 
between the mesosellar behavior, involving three distinct directions of strain, and 
the mesohelical behavior, involving a plane in which material simultaneously rotates 
and, possibly, uniformly strains. Quantities E and A that we defined simplify 
classification of domain points, side-stepping explicit calculation of eigenvalues of 
the flow. The eigenvectors of Vi/’t also have a role in the foregoing analysis. A 
complex conjugate pair determines a plane whose normal is the vector of finite-time 
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rotation. A real eigenvector indicates the direction of finite-time stretching under 
the map. 

Rotation of the material has a much richer presence in 3D than in 2D classifica¬ 
tion. In two dimensions, any non-uniformity in strain, accompanied by rotation or 
not, manifests itself as a hyperbolic deformation. In other words, only rigid-body 
rotations are highlighted as non-hyperbolic in two dimensions, in addition to shear. 
In three dimensions, rotation in a plane can be accompanied by strain in the normal 
direction: the type of that strain distinguishes between three mesohelical classes: 
[—I- -I-], [-h] and neutral. 

Consequently, the invariant vortices in the steady ABC flow initially show a 
significant amount of hyperbolicity in them, indicating that layers within them 
move at different speeds in the axial direction (Figure 3); however, over the longer 
timescales they appear neutrally mesohelical, which corresponds to interpretation 
of them as rigid rotating structures. This suggests that the separation between 
different layers of vortices is asymptotically sub-exponential. 

In the unsteady ABC flow discussed in Section 5.2.3, material is significantly less 
coherent as the unsteadyness destroys long-term invariant structures. Nevertheless, 
for interval lengths T = 1, 5, while the unsteady term A{t) = A-|-Df sin(f) is of the 
same order of magnitude as the steady A, the distributions of mesohyperbolic areas 
show loose correspondence between steady and unsteady ABC flows. However, as 
T —10, the regions that turned neutrally mesohelical in the steady case are instead 
replaced by a mixture of hyperbolic mesosellar regions, indicative of chaotic mixing. 

We have demonstrated that visualization of mesochronic classes corresponds with 
well-known behavior of the fluid-like ABC flow. In particular, it is interesting to 
see how well the mesohelical regions correspond with the vortex zones in which 
Kolmogorov-Arnold-Moser-type structures are known to exist. At this point, the 
theory is immediately applicable to kinematic analysis of the geometric structure 
in chaotic advection in fluids; however, numerical algorithms used in this paper 
serve the proof-of-concept purpose, and more accurate and efficient methods for 
computation of the flow Jacobian [7] can and should be readily used, if applicable. 

Note that the mesochronic classification is invariant to Galilean transformations 
but not to rotating frame transformations [30] , it therefore discovers transport prop¬ 
erties for three-dimensional incompressible flows which are observed in the given 
frame of reference the system is in. It is interesting to consider the relationship 
between the notion of exponential dichotomy [10] and mesohyperbolicity discussed 
here. Mesohyperbolicity (or -ellipticity, -helicity) are notions that are valid on fi¬ 
nite time intervals. They are best used as tools to study the bifurcations of local 
dynamics of a trajectory in time, when changing time intervals are selected, i.e. 
mesohyperbolicity is a function of two parameters, the beginning and final time. It 
appears plausible that, adapting the techniques used in [53], it can be shown that 
if a trajectory is not mesohyperbolic for any interval [ti,t 2 ] inside the interval of 
interest [Ti,r 2 ], then there is no (appropriately defined [53]) finite-time exponential 
dichotomy on that interval. Thus, the concepts defined here have the potential to 
parametrize finite-time stability properties. 

The true test of any method for detection of geometric structures is its usefulness 
to the applied communities, e.g., physical oceanography, and flow engineers. We 
therefore plan to apply the technique to the unsteady testbed flows in [4] , to transi¬ 
tory systems [49] and to more physically-relevant models in order to further verify 
practical use of the mesochronic classification, beyond confirmations obtained for 
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the 2D case. Furthermore, it remains to be understood if the highlighted quantities, 
e.g., the trace, the cofactor trace, and determinant of the mesochronic Jacobian, 
are also useful in dynamic analysis of turbulent fluid flows. 
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Appendix A. Cubic polynomials and 3x3 matrices 


In this section we review formulas and notation for 3x3 matrices. Let A denote 


a 3 X 3 matrix over the field M, A = 
A is 


ail 0-12 0-13 
021 022 023 
O31 O32 O33 


The characteristic polynomial of 


(24) 


p{X) := det(A • Id — A), 


where Id is the unit matrix. The zeros of p(A) are the eigenvalues of A. Define the 
following quantities: 


:=trA, d^:=detA, 

TOA := tr Cof A = det [ ] 

+ det[““ “-]+det[“- 



MESOCHRONIC CLASSIFICATION IN 3D VECTOR FIELDS 


35 


where Cof A is the cofactor matrix of A, containing signed minors.^ By expanding 
equation (24) we can verify that 

(25) p{X) = — <aA^ + ttiaX — (Ia- 

For a polynomial with real coefficients, the zeros Ai,A 2 ,A 3 are either all real, or 
two of them appear as a complex conjugate pair, in which case we denote them by 
Ar, Ac, Ac- By convention, we will assume that the imaginary part 3Ac > 0. 

The discriminant of a cubic polynomial /(A) = A^ + / 2 A^ + /lA + /o is given by 

-Dif) ■= I 8 / 2 / 1/0 - 4 /I /0 + /l/i 
- 4/3 - 27/J, 

its sign determines whether the zeros of / lie on the real axis or not (see Ref. [22], 
§12.1.B). 

Using the notation from (25), we see that /2 = —Ia, fi = rriA, fo = dA, therefore 
Dip) =18tAmAdA - ^t\dA 

+ — 27d\. 

It then holds (see Ref. [34], §10.3, Ex. 10.14) that: 

• T>{p) > 0, the zeros are Ai, A 2 , A 3 G M, 

• 'D{p) = 0 , the zeros are Ai, A 2 = A 3 G K, 

• 'D{p) < 0, the zeros are A^ G K, Ac, Ac G C. 

The following also holds in general for eigenvalues Ai. 2,3 G C, which can be seen 
by expanding (A — Ai)(A — A 2 )(A — A 3 ) and comparing to (25): 

dA — det A = A 1 A 2 A 3 , 
ruA = tr Cof A = Ai A2 + A1A3 + A2A3, 
tA — tr A = Ai + A2 + A3. 

Appendix B. Differential equation for the mesochronic Jacobian 

If an initial condition x G M and initial time to are fixed, the mesochronic Jaco¬ 
bian matrix V/r(x) satisfies a matrix-valued ODE in the variable r, the length of 
the averaging interval. To evaluate the mesochronic Jacobian for the purposes of 
classifying (to, x) into one of the classes described in Section 3, we numerically solve 
a particular matrix-valued initial value problem and evaluate its solution at r = T. 

To derive the ODE for the mesochronic Jacobian, start from the integral expres¬ 
sion for the time-T map and compute its Jacobian: 

nto+T 

'lpr{x) = X+ f{to + t,ll)t{x))dt 

Jto 

pto + T 

Viprix) =ld+ / 77f{to + t,'tfjtix)) ■'\/ijjt{x)dt 

Jto 

Vl/(.(x) = Vf{to+T,'ljjr{x)) ■ Vtprix), 
where ' denotes d/dr. 

We now substitute relation (9), 77'tp^{x) = ld + tW fr{x), linking the time-r map 
'ipT and the mesochronic Jacobian V/t. To simplify notation, we use M(r) := 


^The transpose of Cof A is the adjugate matrix of A. 
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fr{x) and A{t) := V/(io + T,ipr{x)) for the mesochronic and advected velocity 
field Jacobian, respectively. 

(Id + = A{t) ■ (Id + tM{t)) 

(26) M(t) + tM{t) = A{t) + tA{t) ■ M(t) 

M{t) = {^A(t) — It + A{t) ■ M(t). 

The initial condition for the matrix ODE is set at r = 0 when the mesochronic 
Jacobian Vfr{x) is identical to the velocity field Jacobian V/o(a:) = Wf{to,x), 
therefore 

(27) M(0) = ^(0) and M(0) = A(0)^ 

where the last calculation is obtained from (26) by evaluating the first-order expan¬ 
sion of A{t) — M{t) at t = 0. 

Appendix C. Implementation of the mesochronic classification 

In what follows we provide the basic algorithm which we used to produce a 
numerical implementation of mesochronic classification, used to generate images 
analogous to Figures 7, and 9. The core task is to assign a mesochronic class 
(Definition 2.4), corresponding to the time interval [to,to + T], to a fixed initial 
condition x G 

This task can be split into the following sequence of stages.^ 

(1) Compute the trajectory segment x{t) for t G + T] where x{to) = x. 

(2) Evaluate the Jacobian matrix A{t) := V/(to -|- r, x{to + r)) of the velocity 
field along the trajectory x{t). 

(3) Evaluate the mesochronic Jacobian matrix V/,- at the endpoint r = T by 
integrating the initial value problem given by (26) and (27): 

M(t) = (A(t) — M{t)^It -|- A{t) ■ M{t), 

M(0) = A(0), M(0) = A(0)^ 

(4) Compute the determinant = det M{T), cofactor trace = tr Cof M{T), 
evaluate A and E, (10). 

(5) Based on the signs of A and S assign the mesochronic class to the pair of 

initial condition and time interval (p, + T]) using Theorem 3.2. 

Remark to Step (1). We assume that the velocity field f{t,x) can be evaluated 
on the entire domain of interest and throughout the entire compact time interval 
/ G M. If the velocity field is known only on a grid of points {t, x), then interpolation 
of the velocity field is likely needed. 

Remark to Step (2). We assume the knowledge of the Jacobian matrix of the 
vector field evaluated along the trajectory segment [toj^o + t], termed A{t) := 
^f{to + x,x{to -|- t)). If the velocity field was known analytically, it might be 
possible to express V/ analytically as well and evaluate it along points computed 
in Step 1. If not, V/ can be numerically approximated using a central spatial- 
difference scheme with a spatial step 6. On finite-precision computers, S cannot 
be taken arbitrarily small, due to finite-precision effects. There will always exist 
an optimal, non-zero S, which depends on the machine precision and magnitude of 
higher derivatives of / (see Ref. [50], §8). 

^Dependence on values of x, tg, tg + T is omitted, for shortness. 
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Remark to Step (3). Equation (26) is a linear, matrix-valued ODE where the 
vector field Jacobian A{t), computed in Step 2, comes in as both inhomogeneity 
and as the parametric term. This ODE can be discretized using one of the standard 
time-stepping schemes, on the same time points used for discretization of x{t) in 
Step (1). The examples in this paper were computed using an explicit Adams- 
Bashforth stepping scheme with a fixed time step. 

Since each of the presented steps involves some degree of numerical approxima¬ 
tion, the mesochronic Jacobian M(T) will contain numerical noise. As the deriva¬ 
tion of the mesochronic classification criteria hinges on the assumption of incom¬ 
pressibility of the flow, i.e., (5), the necessary criterion for accuracy of the numeri¬ 
cal approximation is that the numerical compressibility tr M(r) -|- r • tr Cof M(t) + 
det M(r) « 0. If the numerical compressibility is significantly larger than 0 when 
evaluated using numerical M{t), this can be taken as an indication of significant 
numerical errors in the mesochronic Jacobian and, consequently, in mesochronic 
classification. 
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